#include 'debug.inc'      
      
      SUBROUTINE SOIL_MPI &
      & (u,v,qv,pt,pts,phi,phis,pp, &
      & phisuf,phig,psuf, &
      & rg,rat,pl, &
      & tsnoi,t2noi,tsurw, &
      & wgnoi,w2noi,wrnoi, &
      & xlake,watdepth,tslake, &
      & alb,z0w,z0hw, &
      & veg,xlai,rsmin, &
      & emis,z0,z0h, &
      & fcveg, &
      & xice, snowini, &
      & rra,rn,h,le, &
      & gsolo,xdd2, &
      & cdm,tsurf,temmin,temmax, &
      & evapor,zonirr,uvsuf, &
      & deltaz,hsuf, &
      & Srad_surf,Lrad_surf, &
      & SraddirUV_surf,SraddirPAR_surf,SraddirIR_surf, &
      & SraddifUV_surf,SraddifPAR_surf,SraddifIR_surf, &
      & fmeth, &
      & xugeos,xvgeos,ugeos,vgeos, &
      & sigma0, &
      & xrg,xra,xpl, &
      & ts,t2,wg,w2,wr, &
      & dx,dy,xmount,ymount, &
      & timeho,dt,dtl, &
      & rhoasurf,ptop,xlatit,fcor, &
      & iclay,isand,iveg, &
      & nx,nx1,ny,ny1,ns,ns1, &
      & nstep,ndados,nf, &
      & ifsoil,iflake,ifcmodel,radpar, &
      & nlakecall, nsoilcall, &
      & verbose, &
      & li,lf, &
      & imonth,iday,ihour,iminu,iseco, &
      & fnout,fnmap,prt)

      use ALLOC, only : &
      & g, stefan, akapa, r, pi, &
      & row

      use MPI_VARIABLES, only : &
      & nxs2i, nxs2e, nys2i, nys2e, &
      & nxi, nxe, nyi, nye, &
      & nx0i, nx1e, nx2e, ny0i, ny1e, ny2e, &
      & nxsh, nysh, &
      & comm3d, ndim, size_MPI, &
      & prplcxy, parallel, &
      & request_send_2d, request_recv_2d

      use LAKE_DRIVER_PARAMETERS, only : &
      & DEFPARDRIVER, n_select_call

      implicit none

!      include 'defpardriver_interface.h'
      include 'soil.inc'

      ! Input/output variables
      real(8), intent(inout) :: u(nxs2i:nxs2e,nys2i:nys2e,0:ns1,1:3)
      real(8), intent(inout) :: v(nxs2i:nxs2e,nys2i:nys2e,0:ns1,1:3)
      real(8), intent(in) :: qv(nxs2i:nxs2e,nys2i:nys2e,0:ns1,1:3)
      real(8), intent(in) :: pt(nxs2i:nxs2e,nys2i:nys2e,0:ns1,1:3)
      real(8), intent(in) :: pts(nxs2i:nxs2e,nys2i:nys2e,0:ns1)
      real(8), intent(in) :: phi(nxs2i:nxs2e,nys2i:nys2e,0:ns1)
      real(8), intent(in) :: phis(nxs2i:nxs2e,nys2i:nys2e,0:ns1)
      real(8), intent(in) :: pp(nxs2i:nxs2e,nys2i:nys2e,1:4)
      real(8), intent(in) :: phisuf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: phig(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: psuf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: rg(nxs2i:nxs2e,nys2i:nys2e) 
      real(8), intent(inout) :: rat(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: pl(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: tsnoi(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: t2noi(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: tsurw(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: wgnoi(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: w2noi(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: wrnoi(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: xlake(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: watdepth(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: tslake(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: alb(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: z0w(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: z0hw(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: veg(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: xlai(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: rsmin(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: emis(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: z0(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: z0h(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: fcveg(1:nv2,nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: xice(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: snowini(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: rra(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: rn(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: h(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: le(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: gsolo(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: xdd2(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: cdm(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: tsurf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: temmin(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: temmax(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: evapor(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: zonirr(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: uvsuf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: deltaz(nxs2i:nxs2e,nys2i:nys2e) 
      real(8), intent(in) :: hsuf(nxs2i:nxs2e,nys2i:nys2e) 

      real(8), intent(inout) :: Srad_surf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: Lrad_surf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: SraddirUV_surf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: SraddirPAR_surf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: SraddirIR_surf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: SraddifUV_surf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: SraddifPAR_surf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: SraddifIR_surf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(out) :: fmeth(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(inout) :: xugeos(1:ns,1:5)
      real(8), intent(inout) :: xvgeos(1:ns,1:5)
      real(8), intent(inout) :: ugeos(1:ns)
      real(8), intent(inout) :: vgeos(1:ns)
      real(8), intent(in) :: sigma0(0:ns1)
      real(8), intent(inout) :: xrg(1:ndados) 
      real(8), intent(inout) :: xra(1:ndados)
      real(8), intent(inout) :: xpl(1:ndados)
      real(8), intent(inout) :: ts(1:2), t2(1:2)
      real(8), intent(inout) :: wg(1:2), w2(1:2)
      real(8), intent(inout) :: wr(1:2)

      real(8), intent(in) :: dx, dy
      real(8), intent(in) :: xmount, ymount 
      real(8), intent(in) :: timeho
      real(8), intent(in) :: dt, dtl
      real(8), intent(inout) :: rhoasurf(nxs2i:nxs2e,nys2i:nys2e)
      real(8), intent(in) :: ptop
      real(8), intent(in) :: xlatit
      real(8), intent(in) :: fcor

      integer(4), intent(in) :: iclay(nxs2i:nxs2e,nys2i:nys2e)
      integer(4), intent(in) :: isand(nxs2i:nxs2e,nys2i:nys2e)
      integer(4), intent(in) :: iveg(nxs2i:nxs2e,nys2i:nys2e)

      integer(4), intent(in) :: nx, nx1, ny, ny1, ns, ns1
      integer(4), intent(in) :: nstep 
      integer(4), intent(in) :: ndados 
      integer(4), intent(inout) :: nf
      integer(4), intent(in) :: nlakecall, nsoilcall
      integer(4), intent(in) :: ifsoil, iflake, ifcmodel
      integer(4), intent(in) :: radpar
      integer(4), intent(in) :: verbose
      integer(4), intent(inout) :: li, lf
      integer(4), intent(in) :: imonth, iday, ihour, iminu, iseco

      character(len=80), intent(in) :: fnout, fnmap

      logical, intent(in) :: prt

      ! local variables

      real(8) :: usuf(nxi:nx1e,nyi:ny1e)
      real(8) :: vsuf(nxi:nx1e,nyi:ny1e)
      real(8) :: Sflux(nxi:nx1e,nyi:ny1e)
      real(8) :: wgw(1:2),w2w(1:2),wrw(1:2),tsw(1:2),t2w(1:2)

      real(8), allocatable, save :: hw_2d(:,:)
      real(8), allocatable, save :: xlew_2d(:,:)
      real(8), allocatable, save :: cdmw_2d(:,:)
      real(8), allocatable, save :: tsw_2d(:,:)
      real(8), allocatable, save :: h_2d(:,:)
      real(8), allocatable, save :: le_2d(:,:)
      real(8), allocatable, save :: cdm_2d(:,:)
      real(8), allocatable, save :: ts_2d(:,:)
      real(8), allocatable, save :: fmeth_2d(:,:,:)

      real(8), allocatable, save :: h10_2d(:,:) 
      real(8), allocatable, save :: l10_2d(:,:) 
      real(8), allocatable, save :: ls10_2d(:,:) 
      real(8), allocatable, save :: hs10_2d(:,:) 
      real(8), allocatable, save :: Ts0_2d(:,:)
      real(8), allocatable, save :: Tb0_2d(:,:) 
      real(8), allocatable, save :: Tbb0_2d(:,:) 
      real(8), allocatable, save :: h_ML0_2d(:,:) 
      real(8), allocatable, save :: extwat_2d(:,:)
      real(8), allocatable, save :: extice_2d(:,:) 
      real(8), allocatable, save :: trib_inflow_2d(:,:) 
      real(8), allocatable, save :: Sals0_2d(:,:)
      real(8), allocatable, save :: Salb0_2d(:,:) 
      real(8), allocatable, save :: us0_2d(:,:) 
      real(8), allocatable, save :: vs0_2d(:,:) 
      real(8), allocatable, save :: Tm_2d(:,:) 
      real(8), allocatable, save :: alphax_2d(:,:)
      real(8), allocatable, save :: alphay_2d(:,:) 
      real(8), allocatable, save :: a_veg_2d(:,:) 
      real(8), allocatable, save :: c_veg_2d(:,:) 
      real(8), allocatable, save :: h_veg_2d(:,:)

      data w2w/1.,1./
      data wgw/1.,1./

      real(8) :: vegw = 0.
      real(8) :: xlaiw = 1.
      real(8) :: rsminw = 30.
      real(8) :: emisw = 1.
      real(8) :: rraw = 1.
      real(8) :: xdd2w = 1.

      real(8), save :: pi180
      real(8) :: p, taixiy, qa, tclds, precip
      real(8) :: xmin, xmax, ymin, ymax
      real(8) :: ttt1, ttt2
      real(8) :: rgw, rnw, ratw
      real(8) :: hw, xlew, gsolow, tsfunc, cdmw, lake_area, fetch
      real(8) :: hour
      real(8) :: h10, l10, ls10, hs10
      real(8) :: Ts0, Tb0, Tbb0, h_ML0, Tm
      real(8) :: extwat, extice
      real(8) :: trib_inflow 
      real(8) :: Sals0, Salb0
      real(8) :: us0, vs0
      real(8) :: alphax, alphay
      real(8) :: a_veg, c_veg, h_veg

      real(8), parameter :: sndensfresh = 100. ! Fresh snow density, kg/m**3
      real(8), parameter :: initsnowfall = 0.05/3600. ! Initial snowfall rate, m/s

      real(8), allocatable :: x(:), xx(:), xx1(:)

      real(8), external :: TLAKE
      real(8), external :: TSOIL

      integer(4), allocatable :: xi(:), xxi(:)
      integer(4) :: init_T
      integer(4), save :: iyear = 2000
      integer(4), save :: ifirst = 0
      integer(4) :: isandw = 1 
      integer(4) :: ivegw = 1
      integer(4) :: ix, iy

      character(len=80) :: fext
      character(len=60), save :: outpath
      character(len=10) :: workchar
      character(len=60) :: xchar

      logical, save :: flag_assim = .false.
      logical, save :: flag_print = .false.

#if debug2
      write(workchar,'(i10.10)') nstep
      call ARRGATHERWRI_MPI(qv(nxs2i,nys2i,0,1), &
      & nxs2i,nxs2e,nys2i,nys2e,0,ns1, &
      & 0,nx1,0,ny1,0,ns1,nx1,ny1,ns1, &
      & 'qv1_soil_1'//workchar)
#endif

      if (parallel .and. nstep == 1) then
        prplcxy(1) = 1
        call ARR2DEXCH(u(nxs2i,nys2i,ns-1,3), &
        & nxs2i,nxs2e,nys2i,nys2e, &
        & nx0i,nx2e,ny0i,ny2e, &
        & nxsh,nysh, &
        & comm3d,size_MPI,ndim,prplcxy, &
        & request_send_2d, request_recv_2d)
        prplcxy(3) = 1
        call ARR2DEXCH(v(nxs2i,nys2i,ns-1,3), &
        & nxs2i,nxs2e,nys2i,nys2e, &
        & nx0i,nx2e,ny0i,ny2e, &
        & nxsh,nysh, &
        & comm3d,size_MPI,ndim,prplcxy, &
        & request_send_2d, request_recv_2d) 
      endif

      if (ifirst == 0) then
        pi180 = 4.*atan(1.)/180.
        allocate (hw_2d(nxi:nx1e,nyi:ny1e)) ; hw_2d = 0.d0
        allocate (xlew_2d(nxi:nx1e,nyi:ny1e)) ; xlew_2d = 0.d0
        allocate (cdmw_2d(nxi:nx1e,nyi:ny1e)) ; cdmw_2d = 0.d0
        allocate (tsw_2d(nxi:nx1e,nyi:ny1e)) ; tsw_2d = 0.d0
        allocate (h_2d(nxi:nx1e,nyi:ny1e)) ; h_2d = 0.d0
        allocate (le_2d(nxi:nx1e,nyi:ny1e)) ; le_2d = 0.d0
        allocate (cdm_2d(nxi:nx1e,nyi:ny1e)) ; cdm_2d = 0.d0
        allocate (ts_2d(nxi:nx1e,nyi:ny1e)) ; ts_2d = 0.d0
        allocate (fmeth_2d(nxi:nx1e,nyi:ny1e,1:2)) ; fmeth_2d = 0.d0 ! 1 - methane flux
                                                                     ! from land
                                                                     ! 2 - methane flux
                                                                     ! from lakes

        if (iflake == 1) then
!         Initializing LAKE model                
          call INIT_LAKE(nxe-nxi+3,nye-nyi+3,fnmap,nlakecall*dt)
!          call DEFPARDRIVER &
!           & (h10_outt=h10, l10_outt=l10, ls10_outt=ls10, hs10_outt=hs10, &
!           &  Ts0_outt=Ts0, Tb0_outt=Tb0, Tbb0_outt=Tbb0, h_ML0_outt=h_ML0, &
!           &  extwat_outt=extwat, extice_outt=extice, &
!           &  trib_inflow_outt=trib_inflow, Sals0_outt=Sals0, Salb0_outt=Salb0, &
!           &  us0_outt=us0, vs0_outt=vs0, Tm_outt=Tm, alphax_outt=alphax, alphay_outt=alphay, &
!           &  a_veg_outt=a_veg, c_veg_outt=c_veg, h_veg_outt=h_veg, &
!           &  init_T_outt=init_T)
          allocate (x(1:12),xx(1:n_select_call),xx1(1:n_select_call))
          allocate (xi(1:27),xxi(1:n_select_call))
          call DEFPARDRIVER &
          & (x(1), x(2), x(3), x(4), x(5), x(6), &
          &  h10, xx(:), l10, ls10, hs10, Ts0, &
          &  Tb0, Tbb0, h_ML0, &
          &  extwat, xx1(:), extice, &
          &  x(7), trib_inflow, Sals0, Salb0, x(8), x(9), &
          &  x(10), us0, vs0, &
          &  Tm, alphax, alphay, a_veg, c_veg, &
          &  h_veg, x(11), &
          &  xi(1), xi(2), xi(3), xi(4), xxi(:), &
          &  xi(5), xi(6), &
          &  xi(7), x(12), xi(8), &
          &  xi(9), xi(10), xi(11), xi(12), &
          &  xi(13), xi(14), xi(15), xi(16), xi(17), &
          &  xi(18), xi(19), xi(20), xi(21), xi(22), &
          &  xi(23), &
          &  init_T, xi(24), xi(25), xi(26), &
          &  xi(27), xchar)
          deallocate (x,xx,xx1)
          deallocate (xi,xxi)

          allocate (h10_2d(nxi:nx1e,nyi:ny1e) ) ; h10_2d = h10
          allocate (l10_2d(nxi:nx1e,nyi:ny1e) ) ; l10_2d = l10 
          allocate (ls10_2d(nxi:nx1e,nyi:ny1e) ) ; ls10_2d = ls10
          allocate (hs10_2d(nxi:nx1e,nyi:ny1e) ) ; hs10_2d = hs10
          allocate (Ts0_2d(nxi:nx1e,nyi:ny1e) ) ; Ts0_2d = Ts0
          allocate (Tb0_2d(nxi:nx1e,nyi:ny1e) ) ; Tb0_2d = Tb0
          allocate (Tbb0_2d(nxi:nx1e,nyi:ny1e) ) ; Tbb0_2d = Tbb0
          allocate (h_ML0_2d(nxi:nx1e,nyi:ny1e) ) ; h_ML0_2d = h_ML0
          allocate (extwat_2d(nxi:nx1e,nyi:ny1e) ) ; extwat_2d = extwat
          allocate (extice_2d(nxi:nx1e,nyi:ny1e) ) ; extice_2d = extice
          allocate (trib_inflow_2d(nxi:nx1e,nyi:ny1e) )
          trib_inflow_2d = trib_inflow
          allocate (Sals0_2d(nxi:nx1e,nyi:ny1e) ) ; Sals0_2d = Sals0
          allocate (Salb0_2d(nxi:nx1e,nyi:ny1e) ) ; Salb0_2d = Salb0
          allocate (us0_2d(nxi:nx1e,nyi:ny1e) ) ; us0_2d = us0
          allocate (vs0_2d(nxi:nx1e,nyi:ny1e) ) ; vs0_2d = vs0
          allocate (Tm_2d(nxi:nx1e,nyi:ny1e) ) ; Tm_2d = Tm
          allocate (alphax_2d(nxi:nx1e,nyi:ny1e) ) ; alphax_2d = alphax
          allocate (alphay_2d(nxi:nx1e,nyi:ny1e) ) ; alphay_2d = alphay
          allocate (a_veg_2d(nxi:nx1e,nyi:ny1e) ) ; a_veg_2d = a_veg
          allocate (c_veg_2d(nxi:nx1e,nyi:ny1e) ) ; c_veg_2d = c_veg
          allocate (h_veg_2d(nxi:nx1e,nyi:ny1e) ) ; h_veg_2d = h_veg
      
!         Ice, snow thickness and depth initial distrubution on lakes/ocean from map file
          h10_2d(nxi:nx1e,nyi:ny1e) = watdepth(nxi:nx1e,nyi:ny1e)
          l10_2d(nxi:nx1e,nyi:ny1e) = xice(nxi:nx1e,nyi:ny1e)
!          hs10_2d(nxi:nx1e,nyi:ny1e) = snowini(nxi:nx1e,nyi:ny1e)
          do iy = nyi, ny1e
            do ix = nxi, nx1e
              if (l10_2d(ix,iy) <= 0.d0) then
                hs10_2d(ix,iy) = 0.d0
              else
                hs10_2d(ix,iy) = snowini(ix,iy)
              endif
            enddo
          enddo

          outpath = 'results/lake/'
        endif
        if (ifsoil == 4) then
!         Initializiong INM soil model
          call INMSOIL_INIT &
          & (tsnoi,wgnoi,xdd2,snowini,dt)
        endif
        if (verbose >= 1) write(*,*) 'The soil is initialized'
      endif

      hour = ihour + real(iminu)/60. + real(iseco)/3600.
!
! soil balance
!
      li = 3 - li
      lf = 3 - lf

!     do iy=2,ny
!     do ix=2,nx
!        write(99,111) 'fi,fis',phi(ix,iy,ns-1),phis(ix,iy,ns-1),ix,iy
! 111    format(a6,3x,f15.8,3x,f15.8,4x,2i3)
!     enddo
!     enddo

      if (verbose >= 2) write(*,*) 'Call rad'
      if (radpar == 1) then
        call RAD_MPI(rg,rat,pl,xrg,xra,xpl, &
        & nx,ny,ndados,nf,dt,nstep, &
        & xugeos,xvgeos,ugeos,vgeos,ns,zonirr)
      else
        !pl = 0.d0
      endif

!     Start debug      
!      write(*,*) 'xpl = ', xpl
!      write(workchar,'(i10)') nstep
!      call ARRGATHERWRI_MPI &
!      & (pl,nxs2i,nxs2e,nys2i,nys2e,0,0, &
!      & 0,nx1,0,ny1,0,0,nx1,ny1,ns1,'pl'//workchar)
!      call ARRGATHERWRI_MPI &
!      & (zonirr,nxs2i,nxs2e,nys2i,nys2e,0,0, &
!      & 0,nx1,0,ny1,0,0,nx1,ny1,ns1,'zonirr'//workchar)
!      STOP
!     End debug      

      do iy = max(nyi,2), nye !2, ny
      do ix = max(nxi,2), nxe !2, nx
!       deltaz(ix,iy)=(phis(ix,iy,ns-1)+phi(ix,iy,ns-1))/g-hsuf(ix,iy)     estava +phis
        deltaz(ix,iy) = 0.5*(phis(ix,iy,ns-1)-phis(ix,iy,ns) &
        & +phi(ix,iy,ns-1)-phi(ix,iy,ns))/g
!        write(99,*) 'deltazmain=',deltaz(ix,iy)
!        write(99,*)'fi(ix,iy,ns-1)fis',phi(ix,iy,ns-1),phis(ix,iy,ns-1)
        if (deltaz(ix,iy) < 0.1) then
          write(*,*) 'DELTAZ:',ix,iy,deltaz(ix,iy) &
          & ,phis(ix,iy,ns-1),phi(ix,iy,ns-1),hsuf(ix,iy),phisuf(ix,iy) &
          & ,phis(ix,iy,ns),phi(ix,iy,ns) &
          & ,phig(ix,iy)
          write(0,*) 'nh3dFatalError=Deltaz<0.1 Unstable!'
          STOP
        endif
!
!   uvsuf: wind speed at the first model level
!
        uvsuf(ix,iy) = dsqrt((0.5*(u(ix,iy,ns-1,2)+u(ix-1,iy,ns-1,2)))**2 &
        &  +(0.5*(v(ix,iy,ns-1,2)+v(ix,iy-1,ns-1,2)))**2)
     
        usuf(ix,iy) = 0.5*(u(ix,iy,ns-1,2) + u(ix-1,iy,ns-1,2))
        vsuf(ix,iy) = 0.5*(v(ix,iy,ns-1,2) + v(ix,iy-1,ns-1,2))
!
! taixiy: air temperature on the first model level
!
        p = sigma0(ns-1)*pp(ix,iy,2) + ptop
        taixiy = (p/1.e5)**akapa*(pts(ix,iy,ns-1) + pt(ix,iy,ns-1,2))
!        write(*,*) taixiy
!       write(*,'(2i5,3e14.7)') ix,iy,taixiy,p,pts(ix,iy,ns-1)
        if (ifsoil == 1 .or. ifsoil == 4) then
          qa = qv(ix,iy,ns-1,2)
        else
          qa = 1.e-3
        endif
        
! radiation, geostrophic wind and precipitation update
!

        if (ifsoil == 4 .and. radpar /= 3) then
!         INM soil model contains carbon models, that use
!         solar radiation in two wavelength bands. This information
!         is currently only available from Clirad radiation scheme.
          call OPTWARNSTOP('ifsoil == 4 .and. radpar /= 3')
        endif

        if (ifsoil == 1 .or. ifsoil == 4) then
          if (radpar == 2) then
            tclds = 0.d0
            call RAD_PAR2(taixiy,qa,psuf(ix,iy),tclds, &
            & rat(ix,iy),rg(ix,iy))
!         tclds is the cloud fraction in the grid cell: temporarily 0.d0        
          elseif (radpar == 3) then
            rg(ix,iy)  = Srad_surf(ix,iy)
            rat(ix,iy) = Lrad_surf(ix,iy)
!     &    + emis(ix,iy)*stefan*tsurf(ix,iy)**4                                  VS,06.2007
          endif
          if(verbose >= 2) write(*,*) 'Rad done'
        endif

        ts(li) = tsnoi(ix,iy)
        tsw(li) = tsurw(ix,iy)
        t2(li) = t2noi(ix,iy)
        wg(li) = wgnoi(ix,iy)
        w2(li) = w2noi(ix,iy)
        wr(li) = wrnoi(ix,iy)


        if (ifirst == 0) then
          ts(lf) = tsnoi(ix,iy)
          tsw(lf) = tsurw(ix,iy)
          t2(lf) = t2noi(ix,iy)
          wg(lf) = wgnoi(ix,iy)
          w2(lf) = w2noi(ix,iy)
          wr(lf) = wrnoi(ix,iy)
        endif

        if (xlake(ix,iy) > 0.) then
         if (iflake == 0) then
           tslake(ix,iy) = tlake(ix,iy,timeho)
!          write(*,*) 'Call solonoi_0',ix,iy
           call SOLONOI (-2,isandw,ivegw,vegw &
           & ,xlaiw,rsminw,alb(ix,iy) &
           & ,emisw,z0w(ix,iy),z0hw(ix,iy),rraw &
           & ,taixiy,qa,uvsuf(ix,iy),rgw,ratw,psuf(ix,iy) &
           & ,pl(ix,iy),deltaz(ix,iy) &
           & ,tsw,t2w,wgw,w2w,wrw &
           & ,rnw,hw,xlew,gsolow,dtl,li,lf &
           & ,xdd2w,cdmw,rhoasurf(ix,iy),tslake(ix,iy))
           fmeth_2d(ix,iy,2) = 0. ! ISBA does not calculate methane flux
         elseif (iflake == 1) then
          if (mod(nstep,nlakecall) == 0 .or. nstep == 1) then
!            if (ix == nx1/2 .and. iy == ny1/2) then
!              flag_print = .true.
!            else
              flag_print = .false.
!            endif
!            call LAKE(taixiy,qa,psuf(ix,iy),usuf(ix,iy),vsuf(ix,iy), &
!            & rat(ix,iy),rg(ix,iy),pl(ix,iy),Sflux(ix,iy),deltaz(ix,iy), &
!            & tsw_2d(ix,iy),hw_2d(ix,iy),xlew_2d(ix,iy),cdmw_2d(ix,iy), &
!            & nlakecall*dt,ix-nxi+1,iy-nyi+1,nxe-nxi+3,nye-nyi+3, &
!            & iyear,imonth,iday,ihour,flag_assim,flag_print,outpath)

            lake_area = dx*dy*xlake(ix,iy) ! The area of all lakes in a cell
            fetch = dsqrt(lake_area/pi) ! Fetch is calculated as a
                                        ! radius of a circle which area is that of all lakes in a cell
            Sflux(ix,iy) = 0.
            call LAKE(taixiy,qa,psuf(ix,iy),usuf(ix,iy),vsuf(ix,iy), &
            & rat(ix,iy),rg(ix,iy),pl(ix,iy),Sflux(ix,iy),deltaz(ix,iy), &
            & nlakecall*dt, &
            & h10_2d(ix,iy), l10_2d(ix,iy), ls10_2d(ix,iy), &
            & hs10_2d(ix,iy), Ts0_2d(ix,iy), &
            & Tb0_2d(ix,iy), Tbb0_2d(ix,iy), h_ML0_2d(ix,iy), &
            & extwat_2d(ix,iy), &
            & extice_2d(ix,iy), fcor, trib_inflow_2d(ix,iy), &
            & Sals0_2d(ix,iy), &
            & Salb0_2d(ix,iy), fetch, xlatit, 0.d0, &
            & us0_2d(ix,iy), vs0_2d(ix,iy), Tm_2d(ix,iy), &
            & alphax_2d(ix,iy), &
            & alphay_2d(ix,iy), a_veg_2d(ix,iy), c_veg_2d(ix,iy), &
            & h_veg_2d(ix,iy), lake_area, &
            & tsw_2d(ix,iy),hw_2d(ix,iy),xlew_2d(ix,iy), &
            & cdmw_2d(ix,iy), fmeth_2d(ix,iy,2), &
            & ix-nxi+1,iy-nyi+1,nxe-nxi+3,nye-nyi+3, &
            & iyear,imonth,iday,hour,init_T,flag_assim, &
            & flag_print,outpath)

            rhoasurf(ix,iy) = psuf(ix,iy)/(r*taixiy)
          endif
          hw =      hw_2d(ix,iy)
          xlew =    xlew_2d(ix,iy)
          cdmw =    cdmw_2d(ix,iy)
          tsw(lf) = tsw_2d(ix,iy)
         else
           print*, 'Invalid lake model identifier!'
           stop
         endif
        else
          hw = 0.
          xlew = 0.
          cdmw = 0.
          tsw(lf) = 0.
          fmeth_2d(ix,iy,2) = 0. 
        endif

!       write(*,*) ix,iy
!     & ,tsw,t2w,wgw
!     &    ,rnw,hw,xlew
!     &    ,cdmw,tslake(ix,iy)
        
        if (xlake(ix,iy) < 1.) then
          
          if (ifsoil == 1) then
!           ISBA soil model
            if(iclay(ix,iy) == -1) then
              tsfunc = tsoil(ix,iy,timeho)
            endif
!           write(*,*) 'Call solonoi_1',ix,iy,iveg(ix,iy)
            call SOLONOI(iclay(ix,iy),isand(ix,iy), &
            & iveg(ix,iy),veg(ix,iy) &
            & ,xlai(ix,iy),rsmin(ix,iy),alb(ix,iy) &
            & ,emis(ix,iy),z0(ix,iy),z0h(ix,iy),rra(ix,iy) &
            & ,taixiy,qa,uvsuf(ix,iy),rg(ix,iy),rat(ix,iy),psuf(ix,iy) &
            & ,pl(ix,iy),deltaz(ix,iy) &
            & ,ts,t2,wg,w2,wr &
            & ,rn(ix,iy),h(ix,iy),le(ix,iy),gsolo(ix,iy),dtl,li,lf &
            & ,xdd2(ix,iy),cdm(ix,iy),rhoasurf(ix,iy),tsfunc)
            fmeth_2d(ix,iy,1) = 0. ! ISBA does not calculate methane flux
          elseif (ifsoil == 4) then
            if (mod(nstep,nsoilcall) == 0 .or. nstep == 1) then
!             INM soil model
!              if (nstep > 10*nsoilcall .and. snowini(ix,iy) > 0) then
!               This procedure works instead of specifying
!               snow cover at initial time moment
!                precip = pl(ix,iy) + initsnowfall
!                snowini(ix,iy) = snowini(ix,iy) - &
!                & initsnowfall*row/sndensfresh*nsoilcall*dt
!              endif
              call INMSOIL_NH3D_INTERFACE &
              & (ix,iy,imonth,iday,ihour,iminu,iseco,ifcmodel,xlatit, &
              & rg(ix,iy),rat(ix,iy), &
              & SraddirUV_surf(ix,iy),SraddirPAR_surf(ix,iy),SraddirIR_surf(ix,iy), &
              & SraddifUV_surf(ix,iy),SraddifPAR_surf(ix,iy),SraddifIR_surf(ix,iy), &
              & veg(ix,iy),alb(ix,iy),xdd2(ix,iy),xlake(ix,iy), &
              & taixiy,qa,psuf(ix,iy),pl(ix,iy),usuf(ix,iy),vsuf(ix,iy), &
              & deltaz(ix,iy),nsoilcall*dt, &
              & hsuf(ix,iy),z0(ix,iy), &
              & fcveg(1,ix,iy), &
              & h_2d(ix,iy),le_2d(ix,iy),cdm_2d(ix,iy), &
              & ts_2d(ix,iy),fmeth_2d(ix,iy,1),rhoasurf(ix,iy))
            endif
            h(ix,iy) = h_2d(ix,iy)
            le(ix,iy) = le_2d(ix,iy)
            cdm(ix,iy) = cdm_2d(ix,iy)
            ts(lf) = ts_2d(ix,iy)
          endif
        else
          h(ix,iy) = 0.
          le(ix,iy) = 0.
          cdm(ix,iy) = 0.
          fmeth_2d(ix,iy,1) = 0.
        endif

        h(ix,iy) = h(ix,iy) !*(1.-xlake(ix,iy)) + hw*xlake(ix,iy)
        le(ix,iy) = le(ix,iy) !*(1.-xlake(ix,iy)) + xlew*xlake(ix,iy)
        cdm(ix,iy) = cdm(ix,iy) !*(1.-xlake(ix,iy)) + cdmw*xlake(ix,iy)
        fmeth(ix,iy) = fmeth_2d(ix,iy,1) !*(1.-xlake(ix,iy)) + &
!        & fmeth_2d(ix,iy,2)*xlake(ix,iy)

!       Debug
!        h(ix,iy) = 0.
!        le(ix,iy) = 0.
!       End debug        

!
! ifsoil=2 means no heat or moisture flux, but there is momentum flux.
!

        if (ifsoil == 2) then
          h(ix,iy) = 0.
          le(ix,iy) = 0.
        elseif (ifsoil == 3) then
          if (ix > nx/2) then
            h(ix,iy) = 100.
          else
            h(ix,iy) = -50.
          endif
          le(ix,iy) = 0.
        else

!  cdm=ro*Cd*u
          tsnoi(ix,iy) = ts(lf)
          tsurw(ix,iy) = tsw(lf)
          tsurf(ix,iy) = (1.-xlake(ix,iy))*ts(lf) + xlake(ix,iy)*tsw(lf)
          t2noi(ix,iy) = t2(lf)
          wgnoi(ix,iy) = wg(lf)
          w2noi(ix,iy) = w2(lf)
          wrnoi(ix,iy) = wr(lf)
          temmin(ix,iy) = min(temmin(ix,iy),tsnoi(ix,iy))
          temmax(ix,iy) = max(temmax(ix,iy),tsnoi(ix,iy))
          evapor(ix,iy) = evapor(ix,iy) + le(ix,iy)/2.5e6*dt
        endif
      enddo
      enddo

      call EXTRAH_MPI &
      & (tsurf, nxs2i, nxs2e, nys2i, nys2e, 0, 0, &
      & 1, nx1, 1, ny1, 0, 0)
      call EXTRAH_MPI &
      & (h, nxs2i, nxs2e, nys2i, nys2e, 0, 0, &
      & 1, nx1, 1, ny1, 0, 0)
      call EXTRAH_MPI &
      & (le, nxs2i, nxs2e, nys2i, nys2e, 0, 0, &
      & 1, nx1, 1, ny1, 0, 0)
      call EXTRAH_MPI &
      & (fmeth, nxs2i, nxs2e, nys2i, nys2e, 0, 0, &
      & 1, nx1, 1, ny1, 0, 0)

!      do i = 1, nx
!      if (ifirst == 0) print*, 'In soil' !, cdmw_2d(i,:)
!	  read*
!	enddo

!	print*, 'ra',rat
!	read*
!      print*, 'The first call of soil finished: STOP'                           VS,15.08.2007
!      STOP                                                                      VS,15.08.2007

      ! This output is not implemented for MPI, must set prt = .false.
      if (ifirst == 0 .and. prt) then
        xmin = (2-1.5)*dx-xmount
        ymin = (2-1.5)*dy-ymount
        xmax = (nx-1.5)*dx-xmount
        ymax = (ny-1.5)*dy-ymount

        call wrigrd(wgnoi,xmin,xmax,ymin,ymax,ttt1,ttt2 &
        &  ,0,nx1,0,ny1,2,nx,2,ny,fext(fnout,'wga.grd'))
        call wrigrd(h,xmin,xmax,ymin,ymax,ttt1,ttt2 &
        &  ,0,nx1,0,ny1,2,nx,2,ny,fext(fnout,'h_a.grd'))
        call wrigrd(le,xmin,xmax,ymin,ymax,ttt1,ttt2 &
        &  ,0,nx1,0,ny1,2,nx,2,ny,fext(fnout,'lea.grd'))
      endif

#if debug1
        write(workchar,'(i10.10)') nstep
        call ARRGATHERWRI_MPI &
        & (tsurf,nxs2i,nxs2e,nys2i,nys2e,0,0, &
        & 0,nx1,0,ny1,0,0,nx1,ny1,ns1,'tsurf_soil'//workchar)
        call ARRGATHERWRI_MPI &
        & (h,nxs2i,nxs2e,nys2i,nys2e,0,0, &
        & 0,nx1,0,ny1,0,0,nx1,ny1,ns1,'h_soil'//workchar)
        call ARRGATHERWRI_MPI &
        & (le,nxs2i,nxs2e,nys2i,nys2e,0,0, &
        & 0,nx1,0,ny1,0,0,nx1,ny1,ns1,'le_soil'//workchar)
        call ARRGATHERWRI_MPI &
        & (cdm,nxs2i,nxs2e,nys2i,nys2e,0,0, &
        & 0,nx1,0,ny1,0,0,nx1,ny1,ns1,'cdm_soil'//workchar)
#endif

#if debug2
      write(workchar,'(i10.10)') nstep
      call ARRGATHERWRI_MPI(qv(nxs2i,nys2i,0,1), &
      & nxs2i,nxs2e,nys2i,nys2e,0,ns1, &
      & 0,nx1,0,ny1,0,ns1,nx1,ny1,ns1, &
      & 'qv1_soil_2'//workchar)
      call ARRGATHERWRI_MPI &
      & (deltaz,nxs2i,nxs2e,nys2i,nys2e,0,0, &
      & 0,nx1,0,ny1,0,0,nx1,ny1,ns1,'deltaz_soil'//workchar)
#endif

      if (ifirst == 0)  ifirst = 1

      return
      END SUBROUTINE SOIL_MPI
